/**
 * @file   Interp.cpp
 * @author xingyifan
 * @date   2021/10/25 10:59:00
 *
 * @brief  member functions
 *
 *
 */

#include "../include/Interp.h"

/**
 * @file   Func.cpp
 * @author xingyifan
 * @date   2021/08/15 16:06:02
 *
 * @brief  polynomial
 *
 *
 */

//计算函数值
double Polynomial::caculate(double x)
{
    double y = 0;
    double x_i = 1;
    int degree = get_degree();
    for (int i = 0; i < degree + 1; i++)
    {
        y += coefficient[i] * x_i;
        x_i *= x;
    }
    return y;
}

//求导
Polynomial Polynomial::derivation()
{
    return derivation(0);
}

Polynomial Polynomial::derivation(u_int n)
{
    vector<double> _coefficient;
    int degree = get_degree();
    for (int i = 1; i < degree + 1; i++)
    {
        _coefficient.push_back(coefficient[i] * i);
    }
    Polynomial poly(_coefficient);
    if (n == 0)
    {
        return poly;
    }
    else
    {
        return poly.derivation(n - 1);
    }
}

//不定积分
Polynomial Polynomial::integration()
{
    vector<double> _coefficient;
    int degree = get_degree();
    _coefficient.push_back(0);
    for (int i = 0; i < degree + 1; i++)
    {
        _coefficient.push_back(coefficient[i] / (i + 1));
    }
    Polynomial func(_coefficient);
    return func;
}

//定积分
double Polynomial::integration(double left, double right)
{
    Polynomial poly(integration());
    double y = poly.caculate(left) - poly.caculate(right);
    return y;
}

Polynomial Polynomial::Chebyshev(u_int n)
{
    if (n == 0)
    {
        return Polynomial({1});
    }
    if (n == 1)
    {
        return Polynomial({0, 1});
    }
    return Polynomial({0, 2}) * Chebyshev(n - 1) - Chebyshev(n - 2);
}

/**
 * @file   Interp.cpp
 * @author xingyifan
 * @date   2021/10/25 21:40:42
 *
 * @brief  NewtonInterp
 *
 *
 */

void NewtonInterp::caculateTable(int start, int end)
{
    for (int i = start; i < end; i++)
    {
        for (int j = 0; j <= i; j++)
        {
            if (j == 0)
            {
                tableOfDividedDiffs_.push_back({interpC_.interpSites[interpC_.points[i]][0]});
            }
            else
            {
                double left = interpC_.points[i - j], right = interpC_.points[i];
                if (left != right)
                {
                    tableOfDividedDiffs_[i].push_back((tableOfDividedDiffs_[i][j - 1] - tableOfDividedDiffs_[i - 1][j - 1]) / (right - left));
                }
                else
                {
                    int factor = 1;
                    for (int k = 1; k <= j; k++)
                    {
                        factor *= k;
                    }
                    tableOfDividedDiffs_[i].push_back(1.0 / factor * interpC_.interpSites[interpC_.points[i]][factor]);
                }
            }
        }
    }
}

void NewtonInterp::caculatePoly()
{
    //由差分表求取插值多项式
    Polynomial poly;
    for (int i = 0; i < (int)tableOfDividedDiffs_.size(); i++)
    {
        Polynomial tmp({tableOfDividedDiffs_[i][i]}); //系数
        //多项式连乘
        for (int j = i - 1; j >= 0; j--)
        {
            tmp *= Polynomial({-1 * interpC_.points[j], 1});
        }
        poly += tmp;
    }
    interPoly_ = poly;
}

void NewtonInterp::overwriteInterp(INTERPC interpC)
{
    for (int i = 0; i < (int)interpC.points.size(); i++)
    {
        interpC_.points.push_back(interpC.points[i]);
        vector<double>::iterator it = interpC.interpSites[interpC.points[i]].begin();
        interpC_.interpSites[interpC.points[i]].push_back(*it);
        interpC.interpSites[interpC.points[i]].erase(it);
    }
    tableOfDividedDiffs_.clear();
    caculateTable(0, interpC_.points.size());
    caculatePoly();
}

void NewtonInterp::incrementInterp(INTERPC interpC)
{
    int size = interpC_.points.size();
    int newSize = interpC.points.size();
    for (int i = 0; i < newSize; i++)
    {
        interpC_.points.push_back(interpC.points[i]);
        interpC_.interpSites[interpC.points[i]] = interpC.interpSites[interpC.points[i]];
    }
    caculateTable(size, size + newSize);
    caculatePoly();
}

double NewtonInterp::Neville_Aitken(double x, INTERPC &interpC)
{
    return Neville_Aitken(x, interpC, 0, interpC.points.size() - 1);
}

double NewtonInterp::Neville_Aitken(double x, INTERPC &interpC, int left, int right)
{
    if (right - left == 0)
    {
        return interpC.interpSites[interpC.points[left]][0];
    }
    if (right - left == 1)
    {
        double poly_i1 = interpC.interpSites[interpC.points[right]][0];
        double poly_i = interpC.interpSites[interpC.points[left]][0];
        return recursion(x, interpC.points[left], interpC.points[right], poly_i1, poly_i);
    }
    //分为三部分
    double poly_left = Neville_Aitken(x, interpC, left, right - 2);
    double poly_mid = Neville_Aitken(x, interpC, left + 1, right - 1);
    double poly_right = Neville_Aitken(x, interpC, left + 2, right);

    poly_left = recursion(x, interpC.points[left], interpC.points[right - 1], poly_mid, poly_left);
    poly_right = recursion(x, interpC.points[left + 1], interpC.points[right], poly_right, poly_mid);

    return recursion(x, interpC.points[left], interpC.points[right], poly_right, poly_left);
}

Polynomial NewtonInterp::Neville_Aitken(INTERPC &interpC)
{
    return Neville_Aitken(interpC, 0, interpC.points.size() - 1);
}

Polynomial NewtonInterp::Neville_Aitken(INTERPC &interpC, int left, int right)
{
    if (right - left == 0)
    {
        return Polynomial(vector<double>({interpC.interpSites[interpC.points[left]]}));
    }
    if (right - left == 1)
    {
        Polynomial poly_i1(vector<double>({interpC.interpSites[interpC.points[right]]}));
        Polynomial poly_i(vector<double>({interpC.interpSites[interpC.points[left]]}));
        return recursionPoly(interpC.points[left], interpC.points[right], poly_i1, poly_i);
    }
    //分为三部分
    Polynomial poly_left = Neville_Aitken(interpC, left, right - 2);
    Polynomial poly_mid = Neville_Aitken(interpC, left + 1, right - 1);
    Polynomial poly_right = Neville_Aitken(interpC, left + 2, right);

    poly_left = recursionPoly(interpC.points[left], interpC.points[right - 1], poly_mid, poly_left);
    poly_right = recursionPoly(interpC.points[left + 1], interpC.points[right], poly_right, poly_mid);

    return recursionPoly(interpC.points[left], interpC.points[right], poly_right, poly_left);
}

double NewtonInterp::recursion(double x, double x_i, double x_i_k1, double poly_i1, double poly_i)
{
    return ((x - x_i) * poly_i1 - (x - x_i_k1) * poly_i) / (x_i_k1 - x_i);
}

Polynomial NewtonInterp::recursionPoly(double x_i, double x_i_k1, Polynomial poly_i1, Polynomial poly_i)
{
    return (Polynomial({-1 * x_i, 1}) * poly_i1 - Polynomial({-1 * x_i_k1, 1}) * poly_i) / (x_i_k1 - x_i);
}

void NewtonInterp::Chebyshev(u_int n, double (*f)(double))
{
    for (int i = 0; i < (int)n; i++)
    {
        double point = cos(1.0 * (2 * i + 1) / n * PI);
        interpC_.points.push_back(point);
        interpC_.interpSites[point] = {f(point)};
    }
    caculateTable(0, interpC_.points.size());
    caculatePoly();
}
